Guaranteeing Spatial Uniformity in Diffusively-Coupled Systems 



S. Yusef Shafi * 



August 22, 2012 



^ . Abstract 

We present a condition that guarantees spatially uniformity in the solution trajectories of a diffusively- 
coupled compartmental ODE model, where each compartment represents a spatial domain of components 
interconnected through diffusion terms with like components in different compartments. Each set of 
like components has its own weighted undirected graph describing the topology of the interconnection 
between compartments. The condition makes use of the Jacobian matrix to describe the dynamics of 
each compartment as well as the Laplacian eigenvalues of each of the graphs. We discuss linear matrix 
inequalities that can be used to verify the condition guaranteeing spatial uniformity, and apply the result 
to a coupled oscillator network. Next we turn to reaction-diffusion PDEs with Neumann boundary 
conditions, and derive an analogous condition guaranteeing spatial uniformity of solutions. The paper 
contributes a relaxed condition to check spatial uniformity that allows individual components to have 
their own specific diffusion terms and interconnection structures. 



1 Introduction 

Diffusively coupled models are crucial to understanding the dynamical behavior of a range of engineering and 
biological systems. Spatially distributed systems whose individual compartments interact with each another 
over weighted undirected graphs are naturally modeled as diffusively-coupled systems of ordinary differential 
equations (ODEs). Coupled ring oscillators constitute a class of voltage-controlled oscillators frequently 
used in applications such as clock recovery circuits and disk-drive read channels PQ P] [3]. Understanding 
synchronization behavior in ring oscillator circuits requires the tools of nonlinear analysis. Turning to the 
biological context, one of the major ideas behind pattern formation in cells and organisms is based on 
diffusion-driven instability [3] [5]. The behavior is observed when higher-order spatial modes in a reaction- 
diffusion partial differential equation (PDE) are destabilized by diffusion, resulting in the growth of spatial 
inhomogeneities, and has been studied extensively [BJ [7J [5] [5]. 

We begin our discussion in Section 2 by studying compartmental ODE models, where each compartment 
represents a well-mixed spatial domain wherein like components in different compartments are coupled by 
diffusion [TO]. In contrast to previous work [TJJ where the interconnections between compartments were 
identical for all components, we allow diffusion terms to be unique to each component. In particular, the 
interconnections within each set of like components may be described by its own weighted undirected graph. 
We derive Lyapunov inequality conditions that guarantee spatial homogeneity in the trajectories of each 
component. We then derive convex linear matrix inequality |12j tests as in that can be used to certify 
the Lyapunov inequality conditions. In Section 3, we apply the LMI tests to study the behavior of a coupled 
ring oscillator circuit. We emphasize that each node in a component may have its own set of neighboring 
nodes to which it is diffusively coupled independent of the set of neighbors of other nodes in the same 
compartment. 
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We next turn to reaction-diffusion PDEs with Neumann boundary condition in Section 4, and establish a 
spatial homogeneity condition analogous to the result for coupled oscillators. The condition we derive allows 
for elliptic operators that model diffusion that may vary spatially and between species [13] , whereas previous 
work has studied the special case of the Laplacian operator [2] [TS] [S] [IZ] [H] ■ We conclude in Section 5. 



2 Spatial uniformity in diffusively coupled systems of ODEs 

We begin by considering a compartmental ODE model where each compartment represents a spatial domain 
interconnected with the other compartments over an undirected graph: 

ii.k = f(xi)k + 51 Wij ( X j,k-Xi,k), i = l,...,N. (1) 

The vector Xi £ M" is the state of the i-th compartment, the vector held /(a;,)^ is the fc-th component of the 
vector field f(xi) acting on Xi, the set Afik consists of the neighbors of the A:-th component of compartment i, 

(k) (k) ' 

and the scalar = £ K is a weighting factor. We aggregate the dynamics of each compartment using 
the stacked vector X = [xj . . . xJj] T and represent the interconnections between like components in different 
compartments by a generalized symmetric positive semidefinite graph Laplacian matrix £ M. NxN : 

X = F(X) - r£ L k ® E,)j X, (2) 

where E]~ = e^e^ £ W nxn is the product of the k-ih standard basis vector et multiplied by its transpose, 
and F(X) = [f(xi) T . . . f(xjy) T ] T . In the event that the fc-th set of like components are not interconnected 
with one another, we set Lk = 0. Define A^' as the second smallest eigenvalue of Lk, and note that since 
Lk^N = 0, 

z T L k z > xf ) z T z (3) 
for all z £ R™ with z _L ljv- Let J(x) = \ x denote the Jacobian of f{x) at x. 

Proposition 2.1 Consider the system 0]. Suppose there exists a convex set X £ K n , a positive definite 
matrix P , and a constant e > such that the following conditions hold: 

P - E X 2 ] Ek^j + (j{x) - J2 P±-eI Vx£X (4) 

PE k + E k P h for each k £ {1, . . . , n} with L k ^ 0. (5) 

Then for any pair (i, j) £ {1, . . . , N} x {1, . . . , N} and any index k £ {1, . . . , n}, we have: x^k{t) — — > 
exponentially as t —> oo. □ 

Proof First recall that z T LkZ > X^ z T z for all z _L l^r, and that z T {Lk®I n )z > X^ z T z for all z _L lN®I n . 
Define the following terms: 



1 N 1 



N 

i=l 



X = l N ®x (6) 
X = X -X. 

Since 53j=i x i = 0> ^ holds that X t (1n ® M) = for all matrices M with n rows. The dynamics of X are 
given by: 

X = F(X) - X - LX 

V , (7) 
= F(X) - X - LX, 
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where L = J2k=i Lk® -Efc- We differentiate the candidate Lyapunov function V — \X T [In ® P)X: 

V = X T {I N ® P)(F(X) -X)- X T (I N eg) P)LX 

n ro\ 

= X T (I N ® P)(F(X) -X)- X T J2(Lk ® PP fc )^. 

fe=i 

We observe that 

(L fc ® PP fe ) + (L fc ® PP fc ) T = L fe ® (PP fc + P fc P), (9) 
and that because condition §5§ holds, there exists a matrix Qk such that Q^Qk = \{PEk + E^P). Then 

X T (L k ® PP fe )X = X T (Pv ® Q£)(L fc ® ^n)(/iV <8> Qfc)* Q 

= Vk{ L k® In)Vk, 



where y/s = (Izv ® Qk)X. Because of the orthogonality relation yu -L In ® In and condition <j3j> , it follows 
that 



X T (I N ® P)LX = E y!(L k ® J„)y fc 
fc=l 



fe=l i=l 

'•T 



Recalling that JT t (1jv <8> M) = 0, we take M = P(f(x) — x) and apply the mean value theorem: 



Because condition (j4]) holds, we have: 



(11) 



fc=l 

n 

fe=l 
n 

= EAf ) l T (/ 7V ®PP fe )l 
fc=i 

n iV 

fe=i »=i 

Defining F(X) — In ® /(#) and adding and subtracting X(In ® P)F(X), we have: 

n AT 

V" < X T (/ A r ® P)(P(X) - X) - E E *! PE A 

k=l i=l 

= X T (I N ® P)(P(X) - + X T (Pv ® P)(1jv ® (/(») - *)) 

n AT 

-E A 2 fc) E^ p ^ ( 12 ) 

fe=i s=i 

= X T (/jv ® P)(P(X) - P(X)) + X T (ljv ® P(/(x) - x)) 

n AT 

-E A 2 fe) E^ Pi? ^- 



(AT \ / n N \ 

E C p - /(«)) - E A ^ fc) E 
i=l / \fc=l i=l / 

N / n \ 

= X>f - /(*)) - P E A 2 fc) ^ ^ ( 13 ) 

»=i V k=i / 

= J2 [ xfP \J(x + SXl ) - E A^ fc) P fe ) ds. 

i=i Jo V k=i J 



< i4 > 
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which concludes the proof. 



□ 



The additional condition that the product PE k be symmetric for each k G {1, . . . , n} with L k ^ allows us 
to generalize Proposition ^. 11 as in to handle non-symmetric generalized graph Laplacians (e.g., where 



,(*) 



^ w$). In place of ©, take 



(L k ® PE k ) + (L fc ® P£ fc ) T = (L fc + Ljf) ® (PE k ), 



(15) 



which holds when is symmetric, and define X 2 as the smallest positive number such that Q holds. 

In order to check the conditions of Proposition 12. 1[ we note two corollaries that follow from ([H], Theorems 
2 and 3), where the Jacobian matrix over the convex set X is itself parametrized by a convex set. 



Corollary 2.2 If there exist constant matrices Z\, . . . , Z q and Si, ... , S m such that 

J{x) G conv{Z\. . . . , Z q } + cone{Si, . . . , S m } \t x E X , 
then the existence of a symmetric matrix P satisfying 



(16) 



P Z k -J2\ { 2 k) E k ) + (z k -J2x ( 2 k) E k I P -< 0, k = 1,. 



fe=i 



fc=i 



PSk + S^P r< 0, fc = l,...,m 



i/ie converse is true. 



(17) 



implies condition /or some e > 0. // J(x) is surjective onto conv{Z\, . . . , Z 9 } + cone{Si, . . . , SVn}, t/ien 



□ 



Next, we define a convex box as: 

box{M , Mi, ... , M p } = {M + lj 1 M 1 + ...+ uj p M p \ uj k G [0, 1] for each k = l,...,p}. (18) 



Corollary 2.3 Suppose that J(x) is contained in a convex box: 

J(x) G box{A Q ,Ai,...,Ai} Vi G X, 



(19) 



where A\, . . . , Ai are rank-one matrices that can be written as A; = BiCj ' , with Bi, Ci G K™. // there exists 
a positive definite matrix V with: 



V = 



P ... 

qi 

... qi 



P G R nxn , ft €R, i = 



(20) 



satisfying: 



A - ELl X 2 k)E k B 



C T —I n 



V < 0, 



(21) 



with B — \B\ . . . Bi] and C = [C\ . . . Ci], then the upper left (positive definite) principal submatrix P satisfies 
condition Q) for some e > 0. If I = 1 and the image of X is surjective onto box{Ao, A\}, then the converse 
is true. □ 
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Define the matrices 



An = 



-m 










' 





















-m 





A x = 











A 2 = 
















—173 

























An = 






a 3 /3 3 



Then it follows that J(x) is contained in a convex box: 

J(x) e box{A , Ax, A 2 , A 3 }. 



(26) 
(27) 

While the method of Corollary 12 . 2 1 involves parametrizing a convex box as a convex hull with TP vertices, and 
potentially a prohibitively large linear matrix inequality computation, the problem structure can be exploited 
using Corollary 12. 31 to obtain a simple analytical condition for synchronization of trajectories. In particular, 
the Jacobian of the ring oscillator exhibits a cyclic structure. The matrix for which we seek a V satisfying 
(|2"Tj) is given by: 



M 



A 



J2k=x ^2 
c T 



B 
-I 



X 



(1) 

2 














-OLxfix 





-T]2 - A 2 





<*2^2 














\(3) 





0:3/33 





1 








-1 











1 








-1 











1 








-1 



(28) 

Note that the matrix M exhibits a cyclic structure, and by a suitable permutation G of its rows and columns, 
it can be brought into a cyclic form M = GMG T . Since M is cyclic, it is amenable to an application of the 
secant criterion [TS], which implies that the condition 



-fef afc ^ , < sec 3 (I) 



nf = i(^ + A/) 



holds if and only if M satisfies 



VM + M T V -< 

for some diagonal V >- 0. Pre- and post-multiplying ([30)l by G T and G, respectively, we have 

G T VGM + M T G T VG ~< 0. 



(29) 



(30) 



(31) 



Note that G T VG is diagonal, and so if M is diagonally stable, then M is diagonally stable as well. We 
conclude that if the secant criterion in (I21H is satisfied, then Corollary 12.31 holds, and so Proposition 12.11 
holds, with: 

Xi t k — Xj.k exponentially as t — » 00 (32) 
for any pair (i, j) 6 {1, ... , N} x {1, . . . , N} and any index k 6 {1, 2, 3}. 

We note that the condition for synchrony that we have found recovers Theorem 2 in pQ , which makes use of 
an input-output approach to synchronization |19j . We have derived the condition using Lyapunov functions 
in an entirely different manner from the input-output approach. 



4 Spatial uniformity in reaction-diffusion PDEs 

Consider the connected, bounded domain f2 C W with smooth boundary d£l, spatial variable £ G f2, and 
outward normal vector n(£) for £ G dil. We consider elliptic operators Lfc given by: 



L k u = V • (A fc (0Vu), A k : Q 



(33) 



G 



where the function A k is symmetric and bounded and 3a > such that for all £ S O and for all C = 
(Ci, C2, • ■ • , Cr) S fi, 531 j a ij(£.)0(j — a lC| 2 - We will study the reaction-diffusion equation: 

^ = f(x) + Cx, (34) 

subject to Neumann boundary conditions Vxi(t, £) • n(£) = V£ G dil, where n(£) is a vector normal to dfl, 
x(f,£) € R n , and 

£x=[V-(A 1 (Z)Vx 1 ) ... V • (A n (£)\7x n )} T (35) 

is a vector of elliptic operators with respect to the spatial variable £ applied to each entry of x(t,£). In a 
reaction-diffusion system, x represents a vector of concentrations for the reactants. We do not emphasize 
well-posedness of solutions to reaction-diffusion PDEs: results on existence of solutions to the reaction PDE 
with Ak = dkl for each k can be found in |13j . 

Define Tt{v} = v — v, where 

* = w\L vm - (36) 

Recall the L 2 (f2) inner product 

(u,v) LHn) = f u T (0v(0d£ (37) 



with the norm ||u||i,f(n) = VW ^)l 2 (0)- 

We now recall a result following from the Poincare principle as in [20], which gives a variational characteri- 
zation of the eigenvalues of an elliptic operator. 

Lemma 4.1 Let be the second smallest Neumann eigenvalue of the operator L k as in \33\) defined on 
the connected, bounded domain fi C M. r with smooth boundary dQ and spatial variable £ 6 f2. Let v = v(£,) 
be a function not identically zero in L 2 (Q) with derivatives J^- € L 2 (Q) that satisfies the Neumann boundary 
condition Vi>(£) • n(£) — and satisfies J^vd^ = 0. Then the following inequality holds: 

Vv ■ (A k Vv) d£, > \ { 2 k) [ v 2 d£. (38) 

□ 

We now show that the solutions of (|34|) achieve spatial uniformity under the conditions (J4j)- : 

Proposition 4.2 Consider the system {3$ . Suppose there exists a convex set X C R", positive definite 
matrix P, and constant e > such that the conditions 0)-{5p hold. Then for every classical solution 
x(t, £) : [0, 00) x £1 — > X, I |7r{x(f , £)}| \l 2 (Q) exponentially as t —> 00. □ 



Proof First define x — n{x}. Note that 



^=ir{f(x)}+£x. (39) 



Consider the candidate Lyapunov functional V(x) = h{x, Px)i,2(n)- Differentiating, we have: 

V(x) < (x, Pn{x}) LH a) + (x, P£x) L '{a) ■ (40) 

We consider the expansion: 

n 

(x,PCx) L 2 {n) =^2(x,PE k £x) L 2 {n) , (41) 



fe=i 
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and note that 

(x,PE k £x) L 2 {n) = (x,PE k £ k x) L 2 {n) , (42) 
where the linear operator C k is defined: 

C k x = [L k x k . . . L k x k . . . L k x k ] T . (43) 

From the condition ([5]), we know there exists a matrix Q k such that Q\Q k = \(PE k + E k P). Substituting, 
we have: 

(x,PE k £ k x) L 2( n) = (Q k x,Q k £ k x) L 2( n) = (y k ,£ k y k ) L *(n), (44) 
where y k = Q k x. Consider the following identity: 

V • (y k ,iA k Vy k ,i) = Vy k>i • (A k Vy k ,i) + y M V • (A k Vy k ,i). (45) 

Integrating both sides, noting the Neumann boundary conditions, and applying the divergence theorem, we 
see that the left hand side of the integrated identity is zero. We then have: 

/ Mm v ■ ( A k^y k , t ) d£ = - I Vj/t, • (A fc Vy M ) d£. (46) 
Jn Jn 

Noting that j Q y k e?£ = Q k f n x d£ = 0, we apply Lemma I4TT1 



Vy M • (A fc Vy M ) d£ > A2 / vU d£, (47) 
n Jn 

where is the second Neumann eigenvalue of L k . Substituting, we have: 

n 

(x,PCx) L 2 (n) = y^^y k ,£ k y k )L*(ci) 

k=l 

n 

<-E^fe,ft)p (!! ) (48) 

k=l 
n 

= -^2 X< 2 ' PE kx) L 2 {n) . 
fc=l 

After adding and subtracting f(x) to the first term on the right hand side of f|40[) . we arrive at: 



V < (x, Pf(x) - f(x)) L 2 (n) - *r (2. Pi ^>L 



(£2) 



fc=l , , 

/ N \ (49) 

S,P ( /(!)-/(!) -5^4*^ 



fc=l / / L 2 (fi) 



An application of the mean value theorem to f(x) — f(x) taken together with condition ([4]) gives: 



V < 



J x T p(^J(x + sx)-^\ { 2 k) E k ^xd£ds 



(50) 



< / / --X J £d£d S <-- 7757 V, 



which concludes the proof. 

□ 
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5 Conclusion 



We have derived Lyapunov inequality conditions that guarantee spatial uniformity in the solutions of com- 
partmental ODEs and reaction-diffusion PDEs even when the diffusion terms vary between species. We have 
used convex optimization to develop tests using linear matrix inequalities that imply the inequality conditions, 
and have applied the tests to coupled ring oscillator circuits. In future work, we will study different spatial 
domains with boundary conditions different from the Neumann condition as well as apply the conditions we 
have derived to biological reaction-diffusion networks. 
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